Extrapolating gravitational-wave data from numerical simulations 
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Two complementary techniques are developed for obtaining the asymptotic form of gravitational- 
wave data at large radii from numerical simulations, in the form of easily implemented algorithms. 
It is shown that, without extrapolation, near- field effects produce errors in extracted waveforms 
that can significantly affect LIGO data analysis. The extrapolation techniques are discussed in 
the context of Newman-Penrose data applied to extrapolation of waveforms from an equal-mass, 
nonspinning black-hole binary simulation. The results of the two methods are shown to agree within 
error estimates. The various benefits and deficiencies of the methods are discussed. 
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I. INTRODUCTION 

As numerical simulations of black-hole binaries im- 
prove, the criterion for success moves past the ability 
of a code to merely persist through many orbits of inspi- 
ral, merger, and ringdown. Accuracy becomes the goal, 
as related work in astrophysics and analysis of data from 
gravitational-wave detectors begins to rely more heavily 
on results from numerical relativity. One of the most im- 
portant challenges in the field today is to find and elim- 
inate systematic errors that could pollute results built 
on numerics. Though there are many possible sources of 
such error, one stands out as being particularly easy to 
manage and — as we show — a particularly large effect: the 
error made by extracting gravitational waveforms from a 
simulation at finite radius, and treating these waveforms 
as though they were the asymptotic form of the radiation. 

The desired waveform is the one to which post- 
Newtonian approximations aspire, and the one sought by 
gravitational-wave observatories: the asymptotic wave- 
form. This is the waveform as it is at distances of over 
lO"'^'* M from the system generating the waves. In typical 
numerical simulations, data extraction takes place at a 
distance of order 100 M from the black holes. At this ra- 
dius, the waves are still rapidly changing because of real 
physical effects. Near-field effects [1, 2, 3] are plainly evi- 
dent, scaling with powers of the ratio of the reduced wave- 
length to the radius, {\/r)^} Extraction methods aiming 
to eliminate the infiuence of gauge effects alone (e.g., im- 
proved Regge-Wheeler-Zerilli or quasi-Kinnersley tech- 
niques) will not be able to account for these physical 
changes. 

Even using a rather naive, gauge-dependent extraction 
method, near-field effects dominate the error in extracted 
waves throughout the inspiral for the data presented in 
this paper [2]. For extraction at r = 50 Af, these effects 
can account for a cumulative error of roughly 50% in 



^ We use the standard notation A = A/27r. 



amplitude or a phase difference of more than one radian, 
from beginning to end of a 16-orbit equal-mass binary 
merger. Note that near-field effects should be propor- 
tional to — at leading order — the ratio of A/r in phase 
and (A/r)^ in amplitude, as has been observed previously 
[2, 4]. Crucially, because the wavelength changes most 
rapidly during the merger, the amplitude and phase dif- 
ferences due to near-field effects also change most rapidly 
during merger. This means that coherence is lost be- 
tween the inspiral and merger/ringdown segments of the 
waveform. 

We can see the importance of this decoherence by look- 
ing at its effect on the matched-filtering technique fre- 
quently used to analyze data from gravitational- wave de- 
tectors. Matched filtering [5, 6, 7] compares two signals, 
si(t) and S2{t). It does this by Fourier transforming each 
into the frequency domain, taking the product of the sig- 
nals, weighting each inversely by the noise — which is a 
function of frequency — and integrating over all frequen- 
cies. This match is optimized over the time and phase 
offsets of the input waveforms. For appropriately nor- 
malized waveforms, the result is a number between 
and 1, denoted (si|s2), with representing no match, 
and 1 representing a perfect match. If we take the ex- 
trapolated waveform as Si and the waveform extracted 
at finite radius as S2, we can evaluate the match between 
them. If the extrapolated waveform accurately represents 
the "true" physical waveform, the mismatch (defined as 
1 — (si I S2)) shows us the loss of signal in data analysis if 
we were to use the finite-radius waveforms to search for 
physical waveforms in detector data. 

The waveforms have a simple scaling with the total 
mass of the system, which sets their frequency scale rela- 
tive to the noise present in the detector. In Figs. 1 and 2, 
we show mismatches between finite-radius and extrapo- 
lated data from the Caltech-Cornell group for a range 
of masses of interest to LIGO data analysis, using the 
Initial- and Advanced-LIGO noise curves, respectively, 
to weight the matches. The value of R denotes the coor- 
dinate radius of extraction for the finite-radius waveform. 

The data in these figures is exclusively numerical data 
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FIG. 1: Data-analysis mismatch between finite-radius 
waveforms and the extrapolated waveform for Initial 
LIGO. This plot shows the mismatch between extrapolated 
waveforms and waveforms extracted at several finite radii, 
scaled to various values of the total mass of the binary sys- 
tem, using the Initial-LIGO noise curve. The waveforms are 
shifted in time and phase to find the optimal match. Note 
that the data used here is solely numerical, with no direct 
post-Newtonian contribution. Thus, for masses below 40 Mq, 
this data represents only a portion of the physical waveform. 

from the simulation used throughout this paper, with 
no direct contributions from post-Newtonian (PN) wave- 
forms. However, to reduce "turn-on" artifacts in the 
Fourier transforms, we have simply attached PN wave- 
forms to the earliest parts of the time-domain waveforms, 
performed the Fourier transform, and set to zero all data 
at frequencies for which PN data are used. The match 
integrals are performed over the intersection of the fre- 
quencies present in each waveform, as in Ref. [8]. 

This means that the data used here is not truly com- 
plete for masses below 40 Mq in Initial LIGO and 110 Mq 
in Advanced LIGO, and that a detected signal would ac- 
tually be dominated by data at lower frequencies than 
are present in this data for masses below about 10 Mq. 
These masses are correspondingly larger for shorter wave- 
forms, which begin at higher frequencies. It is important 
to remember that this type of comparison can only show 
that a given waveform (of a given length) is as good as it 
needs to be for a detector. If the waveform does not cover 
the sensitive band of the detector, the detection signal- 
to-noise ratio would presumably improve given a compa- 
rably accurate waveform of greater duration. Thus, the 
bar is raised for longer waveforms, and for lower masses. 

These figures demonstrate that the mismatch can be 
of order 1% when extracting at a radius of i? = 50 M. 
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FIG. 2: Data-analysis mismatch between finite-radius 
waveforms and the extrapolated waveform for Ad- 
vanced LIGO. This plot shows the mismatch between ex- 
trapolated waveforms and waveforms extracted at several fi- 
nite radii, scaled to various values of the total mass of the 
binary system, using the Advanced-LIGO noise curve. The 
waveforms are shifted in time and phase to find the optimal 
match. Note that the data used here is solely numerical, with 
no direct post-Newtonian contribution. Thus, for masses be- 
low 110 Mq , this data represents only a portion of the physical 
waveform. 



For extraction at i? = 225 M, the mismatch is never 
more than about 0.2%. The loss in event rate would 
be — assuming homogeneous distribution of events in 
space — roughly 3 times the mismatch when using a tem- 
plate bank based on imperfect waveforms [9]. Lindblom 
et al. [10] cite a target mismatch of less than 0.5% be- 
tween the physical waveform and a class of model tem- 
plates to be used for detection of events in current LIGO 
detector data.^ Thus, for example, if these numerical 
waveforms were to be used in construction of template 
banks, ^ the waveform extracted at i? = 50M would not 



^ This number of 0.5% results from assumptions about typical 
event magnitude, template bank parameters, and requirements 
on the maximum frequency of missed events. The parameters 
used to arrive at this number are typical for Initial LIGO. 
We emphasize that these waveforms do not cover the sensitive 
band of current detectors, and thus would not likely be used 
to construct template banks without the aid of post-Newtonian 
extensions of the data. Longer templates effectively have more 
stringent accuracy requirements, so the suitability of these ex- 
traction radii would change for waveforms of different lengths. 
In particular, our results are consistent with those of Ref. [8], 
which included non-extrapolated data of shorter duration. 
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be entirely sufficient, in the sense that a template bank 
built on waveforms with this level of inaccuracy would 
lead to an unacceptably high reduction of event rate. The 
waveforms extracted at i? = 100 Ai" and 225 M, on the 
other hand, may be acceptable for Initial LIGO. For the 
loudest signals expected to be seen by Advanced LIGO, 
the required mismatch may be roughly 10~^ [10]. In this 
case, even extraction at i? = 225 M would be insufficient; 
some method must be used to obtain the asymptotic 
waveform. For both Initial and Advanced LIGO, estimat- 
ing the parameters of the waveform — masses and spins of 
the black holes, for instance — requires still greater accu- 
racy. 

Extrapolation of certain quantities has been used for 
some time in numerical relativity. Even papers an- 
nouncing the first successful black-hole binary evolutions 
[11, 12, 13] showed radial extrapolation of scalar physi- 
cal quantities — radiated energy and angular momentum. 
But waveforms reported in the literature have not al- 
ways been extrapolated. For certain purposes, this is 
acceptable — extrapolation simply removes one of many 
errors. If the precision required for a given purpose al- 
lows it, extrapolation is unnecessary. However, for the 
purposes of LIGO data analysis, we see that extrapola- 
tion of the waveform may be very important. 

We can identify three main obstacles to obtaining the 
asymptotic form of gravitational-wave data from numer- 
ical simulations: 

1. Getting the "right" data at any given point, 
independent of gauge effects (e.g., using quasi- 
Kinnersley techniques and improved Regge- 
Wheeler-Zerilli techniques) ; 

2. Removing near-field effects; 

3. Extracting data along a physically relevant path. 

Many groups have attempted to deal with the first of 
these problems.'* While this is, no doubt, an important 
objective, even the best extraction technique to date is 
imperfect at finite radii. Moreover, at finite distances 
from the source, gravitational waves continue to undergo 
real physical changes as they move away from the sys- 
tem [17], which are frequently ignored in the literature. 
Some extraction techniques have been introduced that at- 
tempt to incorporate corrections for these physical near- 
field effects [18, 19, 20]. However, these require assump- 
tions about the form of those corrections, which wc prefer 
not to impose. Finally, even if we have the optimal data 
at each point in our spacetime, it is easy to see that ex- 
traction along an arbitrary (timelike) path through that 
spacetime could produce a nearly arbitrary waveform. 



* See [14] and [15] and references therein for descriptions of quasi- 
Kinnersley and RWZ methods, respectively. Also, an interesting 
discussion of RWZ methods, and the possibility of finding the 
"exact" waveform at finite distances is found in [16]. 



bearing no resemblance to a waveform that could be ob- 
served in a nearly inertial detector. In particular, if our 
extraction point is chosen at a specific coordinate loca- 
tion, gauge effects could make that extraction point cor- 
respond to a physical path which would not represent any 
real detector's motion. It is not clear how to estimate the 
uncertainty this effect would introduce to the waveforms, 
except by removing the effect entirely. 

We propose a simple method using existing data- 
extraction techniques which should be able to overcome 
each of these three obstacles, given certain very basic 
assumptions. The data are to be extracted at a series 
of radii — either on a series of concentric spheres, or at 
various radii along an outgoing null ray. These data 
can then be expressed as functions of extraction radius 
and retarded time using either of two simple methods 
we describe. For each value of retarded time, the wave- 
forms can then be fit to a polynomial in inverse powers 
of the extraction radius. The asymptotic waveform is 
simply the first nonzero term in the polynomial. Though 
this method also incorporates certain assumptions, they 
amount to assuming that the data behave as radially 
propagating waves, and that the metric itself is asymp- 
totically Minkowski in the coordinates chosen for the sim- 
ulation. 

Extrapolation is, by its very nature, a dangerous pro- 
cedure. The final result may be numerically unstable, in 
the sense that it will fail to converge as the order of the 
extrapolating polynomial is increased. This is to be ex- 
pected, as the size of the effects to be removed eventually 
falls below the size of noise in the waveform data. There 
are likely better methods of determining the asymptotic 
form of gravitational waves produced by numerical simu- 
lations. For example, characteristic evolution is a promis- 
ing technique that may become common in the near fu- 
ture [21, 22, 23, 24]. Nonetheless, extrapolation does 
provide a rough and ready technique which can easily be 
implemented by numerical-relativity groups using exist- 
ing frameworks. 

This paper presents a simple method for implement- 
ing the extrapolation of gravitational- wave data from nu- 
merical simulations, and the motivation for doing so. In 
Sec. II, we begin by introducing an extrapolation method 
that uses approximate tortoise coordinates, which is the 
basic method used to extrapolate data in various pa- 
pers [7, 25, 26, 27, 28] by the Caltech-Cornell collabora- 
tion. The method is tested on the inspiral, merger, and 
ringdown waveform data of the equal mass, nonspinning, 
quasicircular 15-orbit binary simulation of the Caltech- 
Cornell collaboration. We present the convergence of the 
wave phase and amplitude as the extrapolation order in- 
creases, and we also compare data extrapolated using 
various extraction radii. In Sec. Ill, we propose a differ- 
ent extrapolation method using the wave phase — similar 
to the method introduced in Ref. [4] — to independently 
check our results, again demonstrating the convergence 
properties of the method. In Sec. IV, we compare the 
extrapolated waveforms of both methods at various ex- 
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trapolation orders, showing that they agree to well within 
the error estimates of the two methods. A brief discus- 
sion of the pitfalls and future of extrapolation is found in 
Sec. V. Finally, we include a brief appendix on techniques 
for filtering noisy data, which is particularly relevant here 
because extrapolation amplifies noise. 

II. EXTRAPOLATION USING APPROXIMATE 
TORTOISE COORDINATES 

There are many types of data that can be extracted 
from a numerical simulation of an isolated source of 
gravitational waves. The two most common methods 
of extracting gravitational waveforms involve using the 
Newman-Penrose ^'4 quantity, or the metric perturba- 
tion h extracted using Regge-Wheeler-Zerilli techniques. 
Even if we focus on a particular type of waveform, the 
data can be extracted at a series of points along the z 
axis, for example, or decomposed into multipole compo- 
nents and extracted on a series of spheres around the 
source. To simplify this introductory discussion of ex- 
trapolation, we ignore the variety of particular types of 
waveform data. Rather, we generalize to some abstract 
quantity /, which encapsulates the quantity to be ex- 
trapolated and behaves roughly as a radially outgoing 
wave. 

We assume that / travels along outgoing null cones, 
which we parametrize by a retarded time t^^t- Along 
each of these null cones, we further assume that / can be 
expressed as a convergent (or at least asymptotic) series 
in 1/1 — where r is some radial coordinate — for all radii 
of interest. That is, we assume 

J./, N f{k){tict) , , 

fe=0 

for some functions /(fe)- The asymptotic behavior of / is 
given by the lowest nonzero f(k)-^ 

Given data for such an / at a set of retarded times, and 
a set of radii {r^}, it is a simple matter to fit the data 
for each value of trot to a polynomial in l/r. That is, for 
each value of trot, we take the set of data {/(trot, fi)} and 
fit it to a finite polynomial so that 

/(tret,r,)c.f:^^#^. (2) 

k=0 

Standard algorithms [29] can be used to accomplish this 
fitting; here we use the least-squares method. Of course, 
because we are truncating the series of Eq. (1) at fc = iV, 
some of the effects from k > N terms will appear at lower 
orders. We will need to choose N appropriately, checking 



^ For example, if / = r\E'4, then /(gj gives the asymptotic behavior; 
if / = 'I'4, then /(j) gives the asymptotic behavior. 



that the extrapolated quantity has converged sufficiently 
with respect to this order. 

A. Radial parameter 

One subtlety to be considered is the choice of r parame- 
ter to be used in the extraction and fitting. For numerical 
simulation of an isolated system, one simple and obvious 
choice is the coordinate radius R used in the simulation. 
Alternatively, if the data is measured on some spheroidal 
surface, it is possible to define an areal radius i?aroai by 
measuring the area of the sphere along with /, and set- 
ting i?aroai = \/ area/47r. Still other choices are certainly 
possible. 

One objective in choosing a particular r parameter is 
to ensure the physical relevance of the final extrapolated 
quantity. If we try to detect the wave, for example, we 
may want to think of the detector as being located at 
some constant value of r. Or, we may want r to asymp- 
totically represent the luminosity distance. These condi- 
tions may be checked by inspecting the asymptotic be- 
havior of the metric components in the given coordinates. 
For example, if the metric components in a coordinate 
system including r asymptotically approach those of the 
standard Minkowski metric, it is not hard to see that an 
inertial detector could follow a path of constant r param- 
eter. 

Suppose we have two different parameters r and f 
which can be related by a series expansion 

r ^ r [1 + a/r + . . .] . (3) 

For the data presented in this paper, we can show that 
the coordinate radius R and areal radius i?aroai are re- 
lated in this way. Introducing the expansion coefficients 
/(J,), we can write 

m.ur)^±^-^^^±^^. (4) 

fc=0 fe=0 

Inserting Eq. (3) into this formula, Taylor expanding, and 
equating terms of equal order fc, shows that /(g) = /(q) 
and /(I) = /(I). Thus, if the asymptotic behavior of / is 
given by /(q) or /(i), the final extrapolated data should 
not depend on whether r or f is used. On the other 
hand, in practice we truncate these series at finite order. 
This means that higher-order terms could "pollute" /(o) 
or f(^iy The second objective in choosing an r parameter, 
then, is to ensure fast convergence of the series in Eq. (2). 
If the extrapolated quantity does not converge quickly as 
the order of the extrapolating polynomial N is increased, 
it may be due to a poor choice of r parameter. 

The coordinate radius used in a simulation may be sub- 
ject to large gauge variations that are physically irrele- 
vant, and hence are not reflected in the wave's behavior. 
That is, the wave may not fall off nicely in inverse powers 
of that coordinate radius. For the data discussed later 
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in this paper, we find that using the coordinate radius 
of extraction spheres is indeed a poor choice, while using 
the areal radius of those extraction spheres improves the 
convergence of the extrapolation. 



B. Retarded-time parameter 

Similar considerations must be made for the choice of 
retarded-time parameter tj-d to be used in extrapolation. 
It may be possible to evolve null geodesies in numerical 
simulations, and use these to define the null curves on 
which data is to be extracted. While this is an interesting 
possibility that deserves investigation, we propose two 
simpler methods here based on an approximate retarded 
time constructed using the coordinates of the numerical 
simulation and the phase of the waves measured in that 
coordinate system. 

Again, we have two criteria for choosing a retarded- 
time parameter. First is the physical suitability in the 
asymptotic limit. For example, we might want the 
asymptotic tict to be (up to an additive term constant 
in time) the proper time along the path of a detector lo- 
cated at constant r. Again, checking the asymptotic be- 
havior of the metric components with respect to tret and r 
should be a sufficient test of the physical relevance of the 
parameters. Second, we wish to have rapid convergence 
of the extrapolation series using the chosen parameter, 
which also needs to be checked. 

As before, we can also show the equivalence of differ- 
ent choices for the tj-ot parameter. Suppose we have two 
different approximations ^ict and tret that can be related 
by a series expansion 



^ret = irct [l + b/r + 



(5) 



Using the new expansion coefficients /(fe), we can write 



/(fe)(^rct) 



fe=0 



E 

fc=0 



/(fc)(^rct) 



(6) 



Now, however, we need to assume that the functions 
/(fc) can be well-approximated by Taylor series. If this 

is true, we can again show that /(q) = /(o) or, if we have 

/(o) — /(o) = 0; that /(j^) = /(x). The condition that / be 
well-approximated by a Taylor series is nontrivial, and 
can help to inform the choice of /. Similarly, the speed 
of convergence of the extrapolation can help to inform 
the choice of a particular tj-ct parameter. While it has 
been shown [.30] that a retarded-time parameter as sim- 
ple as ^rct = T — Ris sufficient for some purposes, we find 
that convergence during and after merger is drastically 
improved when using a somewhat more careful choice. 

Since we will be considering radiation from an iso- 
lated compact source, our basic model for irot comes from 
the Schwarzschild spacetime; we assume that the system 
in question approaches this spacetime at increasing dis- 
tance. In analogy with the time-retardation effect on 



outgoing null rays in a Schwarzschild spacetime [31], we 
define a "tortoise coordinate" by: 
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where Madm is the ADM mass of the initial data.^ In 
standard Schwarzschild coordinates, the appropriate re- 
tarded time would be given by tret = t — r^,. It is not 
hard to see that the exterior derivative dt^ct is null with 
respect to the Schwarzschild metric. 

Taking inspiration from this, we can attempt to ac- 
count for certain differences from a Schwarzschild back- 
ground. Let T and R denote the simulation's coordi- 
nates, and suppose that we extract the metric compo- 
nents g^^ , g^^, and g^^ from the simulation. We seek 
a ^lot {T, R) such that 



dR 



(8) 



is null with respect to these metric components. That is, 
we seek a t^d such that 



dT 



2g 



OR 



('^ ) = . (9) 



We introduce the ansatz t^ct — t — r^, where t is assumed 
to be a slowly varying function of i?,^ and is given by 
Eq. (7) with R in place of r on the right side. If we ignore 
dt/dR and insert our ansatz into Eq. (9), we have 



rprp 



dt_ 

df 



25 



TR 



1 



dt_ 

dfj \1-2Mabm/R 

2 
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RR 



1 



,l-2MADM/i? 

We can solve this for dt/dT: 
1 



0. (10) 



dt 



dT 1 - 2Mabm/R 



gTR^^^gTRY_ gTTgRR 
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TT 



(11) 



® Kocsis and Loeb [32] pointed out that the propagation of a 
roughly spherical gravitational wave should be affected primarily 
by the amount of mass interior to the wave. Because the waves 
from a merging binary can carry off a significant fraction (typi- 
cally a few percent) of the binary's mass, this suggests that we 
should allow the mass in this formula to vary in time, falling by 
perhaps a few percent over the duration of the waveform. How- 
ever, this is a small correction of a small correction; we have 
not found it necessary. Perhaps with more refined methods, this 
additional correction would be relevant. 

More specifically, we need \dt/dR\ <^ \drt/dR\. This condition 
needs to be checked for all radii used, at all times in the simu- 
lation. For the data presented below, we have checked this, and 
shown it to be a valid assumption, at the radii used for extrap- 
olation. 
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Substituting the Schwarzschild metric components shows 
that we should choose the negative sign in the numerator 
of the second factor. Finally, we can integrate (numeri- 
cally) to find 



Now, in the case where g'^^ is small compared to 1, we 
may wish to ignore it, in which case we have 



1 - 2Madm/R 



(13) 



It is not hard to see that this correctly reduces to t = T 
in the Schwarzschild case. 

For the data discussed later in this paper, we make 
further assumptions that = 1 — 2Madm/^, and that 
R = i?aroai- That is, we define the corrected time 




TT 



-1/.9 
1 - 2Madm/R. 



arcal 



dT' (14a) 



and the retarded time 



ret — f^covr 



(14b) 



We find that this corrected time leads to a significant 
improvement over the naive choice of t{T) = T, while no 
improvement results from using Eq. (12). 



C. Application to a binary inspiral 

To begin the extrapolation procedure, we extract the 
(spin-weight s = —2) {I, m) = (2, 2) component of 
^4 data on a set of spheres at constant coordinate 
radius in the simulation.^ In the black-hole binary 
simulations used here (the same as those discussed in 
Refs. [2, 25, 26, 27]), these spheres are located roughly^ 
every AR w lOMi^.^. (where Mi^r is the sum of the ir- 
reducible masses of the black holes in the initial data) 
from an inner radius of i? = TSMin. to an outer ra- 
dius of i? = 225Mi„, where Mirr denotes the total 
apparent-horizon mass of the two holes at the beginning 
of the simulation. This extraction occurs at time steps 
of AT « 0.5Afirr throughout the simulation. We also 



* See Ref. [27] for details of the extraction procedure. We use >I'4 
data here, rather than Regge-Wheeler— ZeriUi data because the 
\E'4 data from this simulation is of higher quality; it appears that 
the RWZ data is more sensitive to changes in gauge conditions 
after the merger. This problem is still under investigation. 

® Explicitly, the extraction spheres are at radii R/M\^^ = 
{75, 85, 100, 110, 120, . . . , 190, 200, 210, 225}, though we find that 
the final result is not sensitive to the exact placement of the 
extraction spheres. 



measure the areal radius, i?aroai, of these spheres by inte- 
grating the induced area element over the sphere to find 
the area, and defining i?aroai = \/ area/47r. This typically 
differs from the coordinate radius R by roughly Mi„/R. 
Because of gauge effects, the areal radius of a coordi- 
nate sphere changes as a function of time, so we measure 
this as a function of time. Finally, we measure the aver- 
age value of g'^"^ as a function of coordinate time on the 
extraction spheres to correct for the dynamic lapse func- 
tion. The areal radius and are then used to compute 
the retarded time tj-et defined in Eq. (14). 

The gravitational- wave data ^4, the areal radius 
^areai, and the lapse N are all measured as functions of 
the code coordinates T and R. We can use these to con- 
struct the retarded time defined in Eq. (14), using i?arcai 
in place of r. This, then, will also be a function of the 
code coordinates. The mapping between (tret, -Rarcai) and 
(T, R) is invertible, so we can rewrite ^'4 as a function of 

and i?arGal- 

As noted in Sec. II B, we need to assume that the ex- 
trapolated functions are well approximated by Taylor se- 
ries. Because the real and imaginary parts of ^'4 are 
rapidly oscillating in the data presented here, we prefer 
to use the same data in smoother form. We define the 
complex amplitude A and phase ((> of the wave: 



i?: 



arcal 



Mi„^4 = ^e" 



(15) 



where A and are functions of tret and i?areai- Note that 
this definition factors out the dominant 1/r behavior of 
the amplitude. This equation defines the phase with an 
ambiguity of multiples of 27r. In practice, we ensure that 
the phase is continuous as a function of time by adding 
suitable multiples of 2-k. The continuous phase is easier 
to work with for practical reasons, and is certainly much 
better approximated by a Taylor series, as required by 
the argument surrounding Eq. (G). 

A slight complication arises in the relative phase offset 
between successive radii. Noise in the early parts of the 
waveform makes the overall phase offset go through mul- 
tiples of 27r essentially randomly. We choose some fairly 
noise-free (retarded) time and ensure that phases corre- 
sponding to successive extraction spheres are matched at 
that time, by simply adding multiples of 27r to the phase 
of the entire waveform — that is, we add a multiple of 27r 
to the phase at all times. 

Extrapolation of the waveform, then, basically consists 
of finding the asymptotic forms of these functions, A and 
as functions of time. We apply the general technique 
discussed above to A and 0. Explicitly, we fit the data to 
polynomials in 1 / RatcaX for each value of retarded time: 



N 



^(^rct J J^aroal) — ^ ^ 



^(fe)(*rot) 
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The asymptotic waveform is fully described by A^q-^ and 
'/'(o) ■ When the order of the approximating polynomials is 
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important, we will denote by A^r and 4>n the asymptotic 
waveforms resulting from approximations using polyno- 
mials of order N . 

We show the results of these extrapolations in the fig- 
ures below. Figs. 3 through 5 show convergence plots 
for extrapolations using orders N — 1-5. The first two 
figures show the relative amplitude and phase difference 
between successive orders of extrapolation, using the cor- 
rected time of Eq. (14). Here, we define 
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and 
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When comparing waveforms extrapolated by polynomi- 
als of different orders, we use Ni, — Na + 1. Note that 
the broad trend is toward convergence, though high- 
frequency noise is more evident as the order increases, 
as we discuss further in the next subsection. The peak 
amphtude of the waves occurs at time iret/Mrr ~ 3954. 
Note that the scale of the horizontal axis changes just be- 
fore this time to better show the merger/ringdown por- 
tion. We see that the extrapolation is no longer conver- 
gent, with differences increasing slightly as the order of 
the extrapolating polynomial is increased. The oscilla- 
tions we see in these convergence plots have a frequency 
equal to the frequency of the waves themselves. Their 
origin is not clear, but may be due to numerics, gauge, 
or other effects that violate our assumptions about the 
outgoing- wave nature of the data. It is also possible that 
there are simply no higher-order effects to be extrapo- 
lated, so low-order extrapolation suffices. 

Figure 5 shows the same data as in Fig. 4, except that 
no correction is used for dynamic lapse. That is, for 
this figure (and only this figure) , we use t^^t = T — r^,, 
where T is simply the coordinate time. This demon- 
strates the need for improved time-retardation methods 
after merger. Note that the extrapolated data during 
the long inspiral is virtually unchanged (note the dif- 
ferent vertical axes). After the merger — occurring at 
roughly t^et/Mi„ — 3954 — there is no convergence when 
no correction is made for dynamic lapse. It is precisely 
the merger and ringdown segment during which extreme 
gauge changes are present in the data used here [27]. 
On the other hand, the fair convergence of the corrected 
waveforms indicates that it is possible to successfully re- 
move these gauge effects. 



D. Choosing the order of extrapolation 

Deciding on an appropriate order of extrapolation to 
be used for a given purpose requires balancing compet- 
ing effects. As we see in Fig. 3, for example, there is 
evidently some benefit to be gained from using higher- 
order extrapolation during the inspiral; there is clearly 
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FIG. 3: Convergence of the amplitude of the extrap- 
olated ^'4, with increasing order of the extrapolat- 
ing polynomial, A^. This figure shows the convergence of 
the relative amplitude of the extrapolated Newman-Penrose 
waveform, as the order A'^ of the extrapolating polynomial is 
increased. (See Eq. (16).) That is, we subtract the ampli- 
tudes of the two waveforms, and normalize at each time by 
the amplitude of the second waveform. We see that increas- 
ing the order tends to amplify the apparent noise during the 
early and late parts of the waveform. Nonetheless, the broad 
(low-frequency) trend is towards convergence. Note that the 
differences decrease as the system nears merger; this is a first 
indication that the extrapolated effects are due to near-field 
influences. Also note that the horizontal axis changes in the 
right part of the figure, which shows the point of merger, and 
the ringdown portion of the waveform. After the merger, the 
extrapolation is nonconvergent, though the differences grow 
slowly with the order of extrapolation. 



some convergence during inspiral for each of the orders 
shown. On the other hand, higher-order methods amplify 
the apparent noise in the waveform. Moreover, late in 
the inspiral, and on into the merger and ringdown, the 
effects being extrapolated may be present only at low or- 
ders; increasing the extrapolation order would be useless 



So-called "junk radiation" is a ubiquitous feature of initial data 
for current numerical simulations of binary black-hole systems. 
It is clearly evident in simulations as large- amplitude, high- 
frequency waves that die out as the simulation progresses. While 
it is astrophysically extraneous, it is nevertheless a correct result 
of evolution from the initial data. Better initial data would, 
presumably, decrease its magnitude. This is the source of what 
looks like noise in the waveforms at early times. It is less appar- 
ent in h data than in \I'4 data because ^4 effectively amplifies 
high-frequency components, because of the relation <E'4 —h. 
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FIG. 4: Convergence of the phase of the extrapolated 
4'4, with increasing order of the extrapolating poly- 
nomial, N. This figure is the same as Fig. 3, except 
that it shows the convergence of phase. Again, increasing 
the extrapolation order tends to amplify the noise during the 
early and late parts of the waveform, though the broad (low- 
frequency) trend is towards convergence. The horizontal-axis 
scale changes just before merger. 



as higher-order terms would simply be fitting to noise. 

The optimal order depends on the accuracy needed, 
and on the size of effects that need to be eliminated from 
the data. For some applications, little accuracy is needed, 
so a low-order extrapolation (or even no extrapolation) 
is preferable. If high-frequency noise is not considered 
a problem, then simple high-order extrapolation should 
suffice. Of course, if both high accuracy and low noise 
are required, data may easily be filtered, mitigating the 
problem of noise amplification. (See the appendix for 
more discussion.) There is some concern that this may 
introduce subtle inaccuracies: filtering is more art than 
science, and it is difficult to establish precise error bars 
for filtered data. 



We note that — as expected from investigations of near-field ef- 
fects [1, 2, 3] — the second-order behavior of the amplitude greatly 
dominates its first-order behavior [4], Thus, there is no improve- 
ment to the accuracy of the amplitude when extrapolating with 
A'' = 1; it would be better to simply use the data from the largest 
extraction radius. 




FIG. 5: Convergence of the phase of ^'4, extrapolated 
with no correction for the dynamic lapse. This figure 
is the same as Fig. 4, except that no correction is made to ac- 
count for the dynamic lapse. (See Eq. (14) and surrounding 
discussion.) Observe that the convergence is very poor after 
merger (at roughly trct/Afiir — 3954). This corresponds to 
the time after which sharp features in the lapse are observed. 
We conclude from this graph and comparison with the previ- 
ous graph that the correction is crucial to convergence of \I'4 
extrapolation through merger and ringdown. 

E. Choosing extraction radii 

Another decision needs to be made regarding the num- 
ber and location of extraction surfaces. Choosing the 
number of surfaces is fairly easy, because there is typ- 
ically little cost in increasing the number of extraction 
radii (especially relative to the cost of — say — running a 
simulation). The only restriction is that the number of 
data points needs to be significantly larger than the order 
of the extrapolating polynomial; more can hardly hurt. 
More careful consideration needs to be given to the loca- 
tion of the extraction surfaces. 

For the extrapolations shown in Figs. 3 and 4, data 
was extracted on spheres spaced by 10 to ISMirr, from 
R = 75Mirr to R = 225Mi„. The outer radius of 225Mi„ 
was chosen simply because this is the largest radius at 
which data exists throughout the simulation; presumably, 
we always want the outermost radii at which the data 
are resolved. In choosing the inner radius, there are two 
competing considerations. 

On one hand, we want the largest spread possible be- 
tween the inner and outer extraction radii to stabilize the 
extrapolation. A very rough rule of thumb seems to be 
that the distance to be extrapolated should be no greater 
than the distance covered by the data. Because the ex- 
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One waveform is extrapolated using only data from small 
radii, R/Mi„ = (50, . . . , 100}. It is clearly not converged, 
and would require higher-order extrapolation if greater 
accuracy is needed. The source of the difference is pre- 
sumably the near-field effect [2], which is proportionally 
larger at small radii. 

Clearly, there is a nontrivial interplay between the 
radii used for extraction and the order of extrapolation. 
Indeed, because of the time-dependence of the various 
elements of these choices, it may be advisable to use 
different radii and orders of extrapolation for different 
time portions of the waveform. The different portions 
could then be joined together using any of various meth- 
ods [7, 33]. 



III. EXTRAPOLATION USING THE PHASE OF 
THE WAVEFORM 



FIG. 6: Comparison of extrapolation of ^'4 using dif- 
ferent sets of extraction radii. This figure compares the 
phase of waveforms extrapolated with various sets of radii. 
All comparisons are with respect to the data set used else- 
where in this paper, which uses extraction radii R/AIin = 
{75, 85, 100, 110, 120, . . . , 200, 210, 225}. The order of the ex- 
trapolating polynomial is TV = 3 in all cases. 

trapolating polynomial is a function of l/R, the distance 
to be extrapolated is l/i?outcr — l/oo = l/i?outor- The 
distance covered by the data is l/i?innor — 1/^outcr, so if 
the rule of thumb is to be satisfied, the inner extraction 
radius should be no more than half of the outer extrac- 
tion radius, -Rinnor ^ -Router/2 (noting, of course, that 
this is a very rough rule of thumb) . 

On the other hand, we would like the inner extraction 
radius to be as far out as possible. Extracting data near 
the violent center of the simulation is a bad idea for many 
reasons. Coordinate ambiguity, tetrad errors, near-field 
effects — all are more severe near the center of the sim- 
ulation. The larger these errors are, the more work the 
extrapolation needs to do. This effectively means that 
higher-order extrapolation is needed if data are extracted 
at small radii. The exact inner radius needed for extrap- 
olation depends on the desired accuracy and, again, the 
portion of the simulation from which the waveform is 
needed. 

We can compare data extrapolated using different sets 
of radii. Figure 6 shows a variety, compared to the data 
used elsewhere in this paper. The extrapolation order is 
N = 3 in all cases. Note that the waveforms labeled 
R/Mi„ = (50, . . . , 100} and R/Mi„ = (100, . . . , 225} 
both satisfy the rule of thumb that the inner radius 
should be at most half of the outer radius, while the 
other two waveforms do not; it appears that violation 
of the rule of thumb leads to greater sensitivity to noise. 



While the tortoise-coordinate method just described 
attempts to compensate for nontrivial gauge perturba- 
tions, it is possible that it does not take account of all 
effects adequately. As an independent check, we discuss 
what is essentially a second — very different — formulation 
of the retarded-time parameter, similar to one first intro- 
duced in Ref. [4]. If waves extrapolated with the two 
different methods agree, then we can be reasonably con- 
fident that unmodeled gauge effects arc not diminishing 
the accuracy of the final result. As we will explain be- 
low, the method in this section cannot be used naively 
with general data (e.g., data on the equatorial plane). In 
particular, we must assume that the data to be extrapo- 
lated consists of a strictly monotonic phase. It is, how- 
ever, frequently possible to employ a simple technique 
to make purely real, oscillating data into complex data 
with strictly monotonic phase, as we describe below. The 
results of this technique agree with those of the tortoise- 
coordinate extrapolation as we show in Sec. IV. 

Instead of extrapolating the wave phase (j) f^nd ampli- 
tude A as functions of time and radius, we extrapolate 
the time tret and the amplitude A as functions of wave 
phase (j) and radius i?aroai- In other words, we measure 
the amplitude and the arrival time to some radius i?aroai 
of a fixed phase point in the waveform. This is the ori- 
gin of the requirement that the data to be extrapolated 
consist of a strictly monotonic phase (/>(irct j ^aroai) (i-e., 
it must be invertible). For the data presented here, the 
presence of radiation in the initial data — junk radiation — 
and numerical noise cause the extracted waveforms to fail 
to satisfy this requirement at early times. In this case, 
the extrapolation is performed separately for each invert- 
ible portion of the data. That is, the data are divided 
into invertible segments, each segment is extrapolated 
separately, and the final products are joined together as 
a single waveform. 



10 



A. Description of the method 

This extrapolation technique consists of extrapolating 
the retarded time and the amphtude as functions of the 
wave phase (j) and the radius i?areai- In other words, when 
extrapolating the waveform, we are estimating the ampli- 
tude and the arrival time of a fixed phase point at infinity. 
Here, we extract the same ^'4, 5^^, and areal-radius data 
used in the previous section. As in the previous method, 
we first shift each waveform in time using ij-et = icorr ~ > 
where tcorr is defined in Eq. (14) and the basic tortoise 
coordinate is defined in Eq. (7) with areal radius as 
the radial parameter. The amplitude and wave phase 
are again defined using Eq. (15), and the phase is made 
continuous as in Sec. II C. Thus, we begin with the same 
data, shifted as with the tortoise-coordinate method. 

Now, however, we change the method, in an attempt 
to allow for unmodeled effects. Instead of extrapolat- 
ing </)(irot, -Rarcai) and A^t^^t, Raroal) , as with the previous 
method, we invert these functions to get iret(^j -Raroai) 
and A{(j), i?areai) as functions of the wave phase (f). In 
other words, we extrapolate the arrival time and the am- 
plitude of a signal to a coordinate radius R for each wave 
phase value. This is done by fitting the retarded time tret 
and the amplitude A data to polynomials in 1 /i?arcai for 
each value of the wave phase: 

A(i?arcal,0)^5]^^^ , (18a) 

fe=0 areal 

i(i?areal,^)^n+^^P^ , (18b) 
^—Q aical 

where the asymptotic waveform is fully described by 
A(o)((/>) and t(o)((?!'). 

With this data in hand, we can produce the asymp- 
totic amplitude and phase as functions of time by plot- 
ting curves in the t-A and t-cj) planes parametrized by 
the phase. In order to be true, single-valued functions, 
we again need monotonicity of the t(o)(^!') data, which 
may be violated by extrapolation. The usable data can 
be obtained simply by removing data from times before 
which this condition holds. 

Choosing the extraction radii and extrapolation order 
for this method follows the same rough recommendations 
described in Sees. II D and II E. Note also that the restric- 
tion that the data have an invcrtible phase as a function 
of time is not insurmountable. For example, data for ^'4 
in the equatorial plane is purely real, hence has a phase 
that simply jumps from to tt discontinuously. However, 
we can define a new quantity 

w{t) = ^i{t) + ii>4{t). (19) 

This is simply an auxiliary quantity used for the extrap- 
olation, with a smoothly varying, invertible phase. The 
imaginary part is discarded after extrapolation. 
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FIG. 7: Convergence of the amplitude of ^'4 extrap- 
olated using the wave phase, with increasing or- 
der of the extrapolating polynomial. This fig- 
ure shows the convergence of the relative amplitude of the 
extrapolated Newman-Penrose waveform extrapolated using 
the wave phase, as the order A'^ of the extrapolating polyno- 
mial is increased. (See Eq. (18).) Increasing the extrapolation 
order tends to amplify the apparent noise during the early and 
late parts of the waveform, but it improves convergence. The 
vertical axis at tret/Mirr ~ 3950 denotes merger. 



B. Results 

In Figs. 7 and 8 we plot the convergence of the rela- 
tive amplitude and phase of the extrapolated (I, m) = 
(2, 2) mode of the ^'4 waveform for extrapolation or- 
ders iV = 1, . . . , 5. A common feature of both plots is 
that during the inspiral, the higher the extrapolation or- 
der, the better the convergence. However, the noise is 
amplified significantly for large orders of extrapolation. 
This method of extrapolation amplifies high-frequency 
noise significantly, compared to the tortoise-coordinate 
method. 

In the inspiral portion, we have a decreasing error in 
the extrapolation of the phase and the amplitude as the 
wavelength of the gravitational waves decreases. In the 
merger/ringdown portion, a more careful choice of the 
radii and order of extrapolation needs to be made. Since 
near-field effects are less significant in the data extracted 
at larger radii, extrapolation at low order {N — 2,3) 
seems sufficient. Data extrapolated at large order {N — 
4, 5) has a larger error in the phase and amplitude af- 
ter merger than data extrapolated at order TV = 2 or 3. 
Moreover, the outermost extraction radius could be re- 
duced, say to i?outcr/-^^iiT = 165 instead of -RoutorZ-^irr = 
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FIG. 8: Convergence of the phase of ^'4 as a function of 
time extrapolated using the wave phase, with increas- 
ing order TV of the extrapolating polynomial. Again, 
increasing the extrapolation order tends to amplify the ap- 
parent noise during the early and late parts of the waveform, 
though convergence is improved significantly. 



225, without having large extrapolation error at late 
times. Using the radius range R/M-^^ — 75, ... , 160 in- 
stead of the range R/Mi„ = 75, ... , 225 would leave the 
extrapolation error during the merger/ringdown almost 
unchanged, while the extrapolation error during the in- 
spiral would increase by about 70%. 

We note that this method allows easy extrapolation of 
various portions of the waveform using different extrac- 
tion radii and orders since — by construction — the wave 
phase is an independent variable. For example, solve 
for the phase value of the merger (/>mcrgcr (defined as the 
phase at which the amplitude is a maximum), then use 
the radius range i?/Mirr = 75, ... , 225 for all phase val- 
ues less than ^merger and the range R/Mirr = 75, . . . , 160 
for all larger phase values. 

This method has been tested also using the coordinate 
radius R and the naive time coordinate T, in place of 
areal radius and corrected time. We found results sim- 
ilar to those discussed in Sec. II. Using the new time 
coordinate icorr instead of the naive time coordinate T 
improved the extrapolation during the merger/ringdown, 
as found in Sec. II. 

As with the previous extrapolation method, increasing 
the extrapolation order gives a faster convergence rate 
of waveform phase and amplitude, but it amplifies noise 
in the extrapolated waveform. To improve convergence 
without increasing the noise, we need a good filtering 
technique for the inspiral data. The junk-radiation noise 



decreases significantly as a function of time, disappearing 
several orbits before merger. However, this noise could 
be reduced by using more extraction radii in the extrapo- 
lation process, or by running the data through a low-pass 
filter. See the appendix for further discussion of filtering. 



IV. COMPARISON OF THE TWO METHODS 

Both methods showed good convergence of the ampli- 
tude and the phase of the waveform as N increased in 
the inspiral portion. (See Figs. 3 and 7 for the ampli- 
tude, and Figs. 4 and 8 for the phase.) The wave-phase 
extrapolation method was more sensitive to noise. In 
the merger/ringdown portion, both methods have similar 
convergence as N increases, especially when the correc- 
tion is taken to account for the dynamic lapse. The use 
of the time parameter tcorr improved the agreement be- 
tween the methods significantly in the merger/ringdown 
portion for all extrapolation orders. Extrapolating at or- 
der = 2 or 3 seems the best choice as the noise and 
phase differences are smallest for these values. 

In Fig. 9, we show the relative amplitude difference be- 
tween data extrapolated at various orders ( A^ = 1 , . . . , 5) . 
There is no additional time or phase offset used in the 
comparison. Ignoring high-frequency components, the 
difference in the relative amplitude is always less than 
0.3% for different extrapolation orders. Even including 
high-frequency components, the differences between the 
two methods are always smaller than the error in each 
method, as judged by convergence plots. In Fig. 10, 
we show the phase difference between the data extrapo- 
lated using both methods. As in the relative amplitude- 
difference plots, the best agreement is achieved during 
the inspiral portion. Ignoring high-frequency compo- 
nents, the difference is less than 0.02 radians for all or- 
ders. In the merger/ringdown portion, the best agree- 
ment between extrapolated waveforms is at order N = 2 
or 3 where the phase difference is less than 0.01 radians. 



V. CONCLUSIONS 

We have demonstrated two simple techniques for 
extrapolating gravitational-wave data from numerical- 
relativity simulations. We took certain basic gauge in- 
formation into account to improve convergence of the ex- 
trapolation during times of particularly dynamic gauge, 
and showed that the two methods agree to within rough 
error estimates. We have determined that the first 
method presented here is less sensitive to noise, and 
more immediately applies to arbitrary wavelike data; this 
method has become the basic standard in use by the 
Caltech-Cornell collaboration. In both cases, there were 
problems with convergence after merger. The source of 
these problems is still unclear, but will be a topic for 
further investigation. 
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FIG. 9: Relative difference in the amplitude of the two 
extrapolation methods as we increase the order of ex- 
trapolation. The best agreement between both methods 
is at low orders of extrapolation, for which the relative dif- 
ference in the amplitude is less than 0.1% during most of the 
evolution. 
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FIG. 10: Phase difference of the two extrapolation 
methods as we increase the order of extrapolation. 

This figure shows the phase difference between waveforms ex- 
trapolated using each of the two methods. The best agree- 
ment between the methods is at orders N = 2 and 3. The 
relative difference in the phase is less than 0.02 radians during 
most of the evolution. 



As with any type of extrapolation, a note of caution is 
in order. It is entirely possible that the "true" function 
being extrapolated bears little resemblance to the ap- 
proximating function we choose, outside of the domain 
on which we have data. We may, however, have reason 
to believe that the true function takes a certain form. 
If the data in question are generated by a homogeneous 
wave equation, for instance, we know that well-behaved 
solutions fall off in powers of l/r. In any case, there is a 
certain element of faith that extrapolation is a reasonable 
thing to do. While that faith may be misplaced, there are 
methods of checking whether or not it is: goodness-of-fit 
statistics, error estimates, and convergence tests. To be 
of greatest use, goodness-of-fit statistics and error esti- 
mates for the output waveform require error estimates 
for the input waveforms. We leave this for future work. 

We still do not know the correct answers to the ques- 
tions numerical relativity considers. We have no ana- 
lytic solutions to deliver the waveform that Einstein's 
equations — solved perfectly — predict would come from a 
black-hole binary merger; or the precise amount of energy 
radiated from any given binary; or the exact kick or spin 
of the final black hole. Without being able to compare 
numerical relativity to exact solutions, we may be leaving 
large systematic errors hidden in plain view. To elimi- 
nate them, we need to use multiple, independent methods 
for our calculations. For example, we might extract ^'4 
directly by calculating the Riemann tensor and contract- 
ing appropriately with our naive coordinate tetrad, and 
extract the metric perturbation using the formalism of 
Regge-Wheeler-Zerilli and Moncrief. By differentiating 
the latter result twice and comparing to ^'4, we could 
verify that details of the extraction methods are not pro- 
ducing systematic errors. (Just such a comparison was 
done in Ref. [28] for waveforms extrapolated using the 
technique in this paper.) Nonetheless, it is possible that 
infrastructure used to find both could be leading to er- 
rors. 

In the same way, simulations need to be performed 
using different gauge conditions, numerical techniques, 
code infrastructures, boundary conditions, and even 
different extrapolation methods. Only when multiple 
schemes arrive at the same result can we be truly con- 
fident in any of them. But to arrive at the same result, 
the waveforms from each scheme need to be processed 
as carefully as possible. We have shown that extrap- 
olation is crucial for highly accurate gravitational wave- 
forms, and for optimized detection of mergers in detector 
data. 
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APPENDIX: FILTERING 

Extrapolating waveforms containing poorly resolved 
high-frequency components amplifies the magnitude of 
the noise in the signal at infinity. One possible solution 
to the problem is to filter out the junk radiation from the 
gravitational waveform. This is possible when the noise 



has higher frequency than the physical data of interest. 
The Matlab function f iltf ilt, using a low-pass But- 
terworth filter with cutoff frequency between the noise 
frequency and the highest gravitational- wave frequency, 
is satisfactory for many uses when applied to the early 
parts of the data. This filtering may be applied to either 
the complex data, or to its amplitude and phase — the 
latter allowing for a lower cutoff frequency. There is also 
a marginal benefit to be gained when the input data are 
filtered before extrapolation, though filtering of the fi- 
nal result is also necessary. It is also possible to fit a 
low-order polynomial to the data, filter the residual, and 
add the filtered data back to the fit; this removes very 
low- frequency components, reducing the impact of Gibbs 
phenomena. 

For the data presented here, we use a sixth-order 
Butterworth filter with a physical cutoff frequency of 
^cutoff = 0.075/Mi„,^^ which is roughly eight times 
the maximum frequency of the physical waveforms ex- 
pected in the filtered region. The filter is applied indi- 
vidually (using the f iltf ilt function) to the amplitude 
and phase data, in turn. Because of remaining Gibbs 
phenomena at late times, we use unfiltered data after 
i,et/Mi„ = 3000. 

One basic diagnostic of the filtering process is to sim- 
ply look at the difference between filtered and unfiltered 
data. If there are low-frequency components in these 
curves, we know the cutoff frequency needs to be raised. 
In Fig. 11, we show the difference in relative amplitude 
(upper panel), and phase (lower panel). Because there 
is no difference between the filtered and unfiltered wave- 
forms on the timescale of the physical gravitational waves 
(> lOOMirr), we conclude that the filter's cutoff fre- 
quency is high enough to retain the physical information. 

On the other hand, to check that the filter's cutoff 
frequency is low enough to achieve its purpose, we look 
at data which previously showed the undesirable high- 
frequency characteristics. In Fig. 12, we show the same 
data as in Fig. 4, when the data are filtered before sub- 
traction. The size of the noise at early times is greatly 
reduced. There are still significant high-frequency fea- 
tures in the plot, though they are much smaller than in 
the unfiltered data. Presumably these features are sim- 
ply so large in the input data that even with the large 
suppression from the filter, they are still noticeable. It 
may be possible to remove them by further decreasing the 
filter's cutoff frequency, though this would require better 
handling of Gibb's phenomena from the beginning and 
end of the wave. We find the present filter sufficient for 
the demonstration purposes of this appendix. 
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